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ABSTRACT 


The differential equations of motion for a coupled two-body, gravity- 
oriented, satellite are derived in Part I. Each of the two bodies is as- 
sumed to be rigid. The satellite configuration is basically an x with a 
damper rod attached so that it is free to move constrained only by a 
spring damper mechanism in a plane fixed with respect to the x plane. 
This leads to a seven degree-of-freedom problem. A FORTRAN II digital 
program is included in Appendix D to provide a numerical solution for 
these equations. 

In Part II one method of orbit adjust and station keeping is studied. 
The satellite parameters are chosen for a medium altitude (approxi- 
mately 6000 nautical miles) orbit. A one pound thrust level and a mis- 
alignment angle of one degree are assumed for the station keeping 
thrusters. The results of the study show that a 0.01 orbit eccentricity 
change can be made in 40 days with tangential pulses of 40 seconds dura- 
tion every 22.4 hours without causing satellite inversion. However, 
tangential pulsing requires sensing the satellite yaw angle and, a signifi- 
cant number of damper rod collisions, with mechanical stops placed at 
one radian on either side of the null position, resulted. These two prob- 
lems are avoided by using radial pulses but, the eccentricity change is 
only 0.0057 after 40 days of correction. 
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THE DEVELOPMENT OF A MATHEMATICAL MODEL AND A 
STUDY OF ONE METHOD OF ORBIT ADJUST AND STATION 
KEEPING GRAVITY-ORIENTED SATELLITES* 


by 

Curtiss C. Barrett 
Goddard Space Flight Center 


INTRODUCTION 

Spacecraft designers are currently faced with the problem of developing, at minimum cost, 
satellites with long lifetimes without sacrificing reliability. Many of these satellites will provide 
services which do not require precise attitude control, such as providing communications links and 
gathering meteorological data. The cost of these satellites can be greatly reduced by having one 
side always facing the earth. This, for example, eliminates the necessity of having directional 
control of television cameras or receiving and transmitting antennas. Gravity-gradient stabiliza- 
tion has this property, and it seems that reliability for long lifetimes is possible, since this is a 
passive technique. 

The concept of gravity- gradient stabilization is not new. Roberson (Reference 1) attributes to 
Laplace the first recognition of this effect and its application in 1785 to the problem of the stability 
of the moon with one face always pointing toward the earth. This tendency of an elongated body to 
align itself with the local gravitational field can be made an actuation control source, if the satel- 
lite configuration is chosen so that the stable attitude in the gravitational field is the desired atti- 
tude. Any system based on this principle would require some method of damping oscillatory 
motion to be practical. Among the early investigators of gravity-gradient stabilization, Paul 
(Reference 2) showed that the librations of an extensible dumbbell satellite could be damped by dis- 
sipating energy through the vibration of an elastic material. About the same time, work was initi- 
ated at the Applied Physics Laboratory of Johns Hopkins University on connecting two masses by 
an elastic spring and energy- dissipating device, a concept which was later incorporated into the 
TRAAC satellite (Reference 3). Kamm (Reference 4) analyzed a satellite to which were hinged, 
with appropriate damping mechanisms, two long rods to provide three-axis damping of librations. 
Tinling and Merrick (Reference 5) showed that similar results could be obtained by a single damper 
rod positioned so as to cause satellite inertial asymmetry with respect to the orbital plane. It was 


•Taken from a dissertation submitted at the Catholic University of America in partial fullfillment of the requirements for the degree of 
Master of Science. 
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shown that the necessary inertial properties could be obtained by erecting fixed rods in the form 
of an X with the payload at the center and a damper rod attached in the horizontal plane. This 
latter concept seems very promising and is currently being considered as a stabilization system 
for a number of iuture satellites. 

The purpose of Part I of this paper is to analyze the motion of the coupled satellite in the 
presence of a central gravitational field. The orientation of the satellite is specified at any time 
by three translational coordinates, three rotational coordinates or attitude deviation angles, and 
an angle which locates the damper rod with respect to the main body of the satellite. These seven 
variables are derived in differential equation form. 

Perturbing forces produced by the non- sphericity of the earth T s gravity field and the presence 
of the sun and moon constantly act to remove the satellite from its desired orbit. Since gravity- 
oriented satellites are expected to endure long lifetimes precisely on station, station keeping be- 
comes a requirement. Thruster misalignment is typical of the many problems to be considered 
and may be even more critical for gravity- oriented satellites because of the large moments of 
inertia involved. 

Part II of this paper presents the results of a study to find a thrust program capable of station 
keeping a gravity-oriented satellite. Disturbing torques considered are those due to thruster mis- 
alignment and those arising naturally from variation of the radius vector in elliptic orbits of small 
eccentricity. 

In Appendix D is a Fortran n digital program which provides the numerical solution to the 
equations derived in Part I. This program was used in the study made in Part II. 


I: DEVELOPMENT OF A MATHEMATICAL MODEL 
FOR A GRAVITY-ORIENTED SATELLITE 


SATELLITE DESCRIPTION 


The gravity- oriented satellite (Figure 1) consists of a central housing from which extend four 
elongated rods in the shape of an X . Coupled to the housing through a spring damping mechanism 
is another rod which is free to rotate in a fixed plane with respect to the plane of the x . Rotation 
is limited by stops on either side of the rest position of the damper rod. In flight, the earth- 
pointing axis of the satellite bisects the acute angles formed by the x, and the plane of the x would 
be skewed to the orbital plane as indicated in the figure by the angle between arms of the x and the 
velocity vector. For the purpose of this report the two individual bodies are assumed rigid. 


The inertial properties of the satellite may 
be varied by changing the angle of the X and by 
varying the angle between the damper rod rota- 
tion plane and the X plane. For the purpose of 
this report these two quantities will be con- 
sidered as constants. 

REFERENCE COORDINATE SYSTEMS 

A complete description of satellite motion 
utilizes several orthogonal right-handed refer- 
ence systems, which are described in this 
section. 

Inertial Coordinate System 

The inertial coordinate system is refer- 
enced to the center of the earth, which is used 
as the origin of the system (Figure 2). The co- 
ordinate axes are: 

X In equatorial plane, along line of equi- 
noxes toward vernal equinox, 

Y In equatorial plane, normal to X, 



Figure 1 -Sketch of gravity-oriented satellite. 
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and 

Z Along earth’s spin axis, toward North 
Pole. 

Geocentric Coordinate System 

The geocentric coordinates, like the iner- 
tial coordinates, are referenced to the center of 
the earth (Figure 3). The coordinate axes are: 

X c Outward along earth- satellite radius 
vector, 

Y c In orbital plane, parallel to satellite 
height circle tangent, 

and 

Z Parallel to satellite meridian circle 

c 

tangent. 


ro NORTH POLE 
z 



POINT OF 
ARIES 

Figure 2— Inertial reference system. 


Attitude Reference System 


The attitude coordinates are referenced to the satellite (Figure 4). The coordinate axes are: 

x 0 Outward along earth- satellite radius vector, 

Y 0 Tangent to satellite orbit path, same sense as vehicle motion, 

and 

Z 0 Normal to satellite orbit plane. 


Body Frame Reference System 

The body frame is referenced to the plane of the X , with the frame center as the origin (Fig- 
ure 5). The coordinate axes are: 

x b In plane of rods, bisecting smaller angle made by rods, 

Y b In plane of rods, normal to X b , 

and 

z b Normal to plane of rods. 

Damper Rod Reference System 

The damper rod is referenced to the rod hinge with the hinge point as origin (Figure 6). The 
coordinate axes are: 
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Y c Tangent to height circle 
Z Tangent to meridian circle 


s i n A cos 4> c 
s in 0 c 


cos A 
0 



Figure 3-Geocentric coordinate system, (a) Geocentric reference frame, shown with satellite origin for clarity; 
rotational properties are not affected, (b) Rotational transformation of coordinates from geocentric to inertial refer- 
ence frame. 


INSTANTANEOUS 




X Q Along radius vector 
Y q In pla ne of orbit 
Z Q Normal to plane of orbit 

Figure 4— Attitude coordinate system, (a) Orbit or attitude reference frame, shown with earth center origin for 
clarity; rotational properties are not affected, (b) Rotational transformation of coordinates from the attitude to the 
inertial reference frame. 
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Figure 5-Body frame coordinate system, (a) Body reference frame, (b) Rotational transformation of coordinates from 

body to attitude reference frame. 



x d Lying along rod length, 

Y d Axis about which rod rotates, 

and 

Z d Defined by T when p, and Pj vanish, 

d~*b 

where the transformation is 



In Figures 4(b), 5(b), and 6, transformation 
matrices accompany the diagrams. At times 
in the report these matrices will be referred 
to merely as ( T ) when they represent an 
obvious transformation of a vector from its 




representation in the [X, Y, z] a system to representation in the [X, Y, z] 0 system, where a and 
/ 8 are coordinate identification subscripts. 


DERIVATION OF EQUATIONS OF MOTION 

Orientation of the satellite is given at any time by the following variables: 

(1) The three translational coordinates, which fix the satellite position with respect to a 
standard inertial coordinate system. 

(2) The three rotational coordinates or attitude deviation angles, designating the satellite 
attitude as measured in an orbital coordinate system or attitude reference system. 

(3) The angle made by the damper rod constrained to rotate in a plane fixed with respect to 
the main body (damper rod angle). 

Translational Equations of Motion 

When viewed in inertial spherical coordinates, the satellite position is determined by the 
parameters r, A, and <f > c . Associated with this triplet are both a velocity vector 


v = [wy] 


and an acceleration vector 



r -i 


_ F r F a 

V 

A = 

1 

o 

< 

< 

< 

< u 

- 

_ m m 

m _ 


where m represents the satellite mass and F r , F A , and F 0 are the external forces acting on the 
system. 


When expressed in spherical coordinates, the components of V and A are 

V r = r , 

V A = r A co s v' , 


v 0 = r^ c . 


(la) 

(lb) 

(lc) 


and 


A r - — = V - rA 2 cos 2 4> q - r4> c 2 , 


(2a) 
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Fa 

a a ' 17 


ft (r 2 Acos 2 <^ c ) , 


* 4 > m 


r i r2 k) + t A2 siri ^c 


Rewriting Equations 2 in terms of the highest derivative results in 


r - (r A 2 cos 2 <fr c + r<£ c 2 ) + ’ 

^ (r 2 Acos 2 <?!> c ) = ^ (rcos^) , 


^ ( r2 ^ c ) 


r A 2 sin 4> cos 4> 


Successive integrations of Equations 3 will lead to the determination of the variables r, A, and <£ c . 
In the present form the equations permit the inclusion of any external force functions known to be 
acting on the system. 

Any second-order ordinary differential equation, when integrated, will require two constants 
of integration. For the definite integrals occurring here, these are the initial conditions, which 
must be determined by other means. The Runge-Kutta method of numerical integration, utilized 
in this report, requires knowledge of the dependent variable and its first derivative at one point 
in order to generate a point-to-point solution for each equation. The total of six initial conditions 
may be obtained (as indicated in Appendix B) by specifying six orbital elements, which may then 
be used to compute the initial conditions. The procedure adopted for computer solution is to 
designate a set of orbital elements, solve for the desired constants, and then integrate the equations 
by quadrature to obtain the time- varying coordinate parameters. Also shown in Appendix B are 
equations necessary for computing new orbital elements from the integrated variables. These or- 
bital changes would be caused by noncentral forces. 


Rotational Equations of Motion 

Once the satellite’s position in space is known, it then remains to determine its attitude with 
respect to a known coordinate system, chosen here to be the attitude reference system shown in 
Figure 4. Such a determination is made easier by the fact that, to an excellent approximation, the 
angular motion of the vehicle about its center of mass can be completely decoupled from its angular 
motion in orbit about the earth. Consequently, the basic equation of attitude motion is simply that 
the time rate of change of angular momentum H is equal to the applied external torque. Consider- 
ing that the body reference system [x b , Y b , z b ] for which the angular momentum is defined is 
itself rotating with respect to the previously noted inertial reference system with angular velocity 
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n = p7 x , Q y , J, the time rate of change of H, as seen by an observer at the origin of the iner- 
tial system, is expressed by 


H c = H + n x H , 

where 



(4) 


(Reference 6). 

For the coupled satellite, Equation 4 may be applied to both the main frame and the damper 
rod, leading to 


d_ 

dt 



yy 

I 


yz 



'fl H 


n h 

*b % 

0 H 

*b y b 



-n 

y b 



for the main frame, and 


and 


i n 

r x . 


M + R 

x d x d 


+ (l r - k) Q Z Q x = K + R y 

y d ' r L/ z d d y d y d 

- (I r - I[)C ft = M + R 


(5) 


(6a) 

(6b) 

(6c) 


for the hinged rod. Note that l r is about the rod axis, and if is normal to the axis. 

In Equations 5 and 6 the term M represents external torques on the configuration, while R rep- 
resents the internal reaction torque(s) between the two sections of the vehicle. It follows that since 
the rod and frame are connected at all times, the reaction torques on the two sections are equal in 
magnitude although opposite in direction. This may be expressed as 



(7) 
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The spring-damper mechanism acts, as mentioned earlier, to damp the rod in its rotation 
about the Y d axis. This restoring torque, denoted by M SD , is equivalent to the reaction torque 
component R y ; if Equations 6a to 6c are rearranged to solve for R x ^, R^ and R z ^ , respectively, 
the results are 




R Mc n i 


and 


- (!,- ' M z 


(8a) 

(8b) 

(8c) 


Were the damper rod to be confined rigidly to the Y d axis, the transformation of angular 
velocity fi from the body coordinate system to the rod coordinate system could be expressed 
simply as fl, - ( b T d ) n t>- However, the motion p 2 of the rod induced relative to the Y d axis by the 
spring coupling must be considered in the transformation. This is seen to lead to the revised 
expression 



(9) 


Formal differentiation of Equation 9 under the rules of matrix algebra yields 



Note that in the actual differentiation of the matrix j the angle p l made by the plane of 
rotation of the rod with the X-frame plane is a constant, and so all d/dt ( p , ) terms vanish. 

At this point it should be explained what is desired: a functional relation between the angular 
velocity vector n = [n , n n 1 and known quantities. This expression, in the form of a 
differential equation or set of equations, will then be numerically integrated to determine the angu- 
lar velocity components. However, these are in turn functionally related to the derivatives of the 
attitude deviation angles a, <£, andi//; thus all that remains is to integrate the resulting set of 
equations. In this manner the attitude deviation angles will have been determined. 

To accomplish this task, Equation 8 is first rewritten in terms of n b and H b ; these latter 
terms come from Equations 9, 10a, and 10c. It should be pointed out that Equation 10b is not 


10 


# 

used at the present time. The revised Equation 8 is substituted for the vector [R Xd , R yj , R z<j J 
in Equation 7. The right-hand term in Equation 7 then replaces |"r x , R , R 1 in Equation 5; 

L b J b b 

the new equation results. 



, I \ 




( a h -n h \ 


u \ 

If , (n b .n b )\ 

XX 

xy xz | 


f x b 


f y b * b % y b | 


Xb 1 


f l \ d b/ 

I 

I I 


n 

+ 

n h -n h 

= 

M 

+ 

f, (a , n. ) 

xy 

yy yz 


, ' 


z b x b x b z b 


y b i 


2 \ b b/ 

V- 

^ 1 



f 

\n h -n h / 

\ x b y b x b / 


VJ 


\^ f 3 (^b’ \)j 


Comparison with Equation 5 reveals that the term (d/dt) [ l ay3 ] is assumed to be of the form 


( 0 0 
0 0 
0 0 



The terms f x , f 2 , and f 3 are written here to show that a functional relationship exists in terms 
of the quantities O b and n b . The actual expansions are not included here, since they are merely 
an intermediate step in the solution of the problem. 


Equation 11 may be rearranged so that all derivatives of the angular velocity appear on one 
side of the equation; this is represented by 



where 


A 


det 


Jxy J x z 

j yy J yz 

Jyz J ZZ 


and the matrix elements K.. = (a..) (cofactor of J Aj ), where i, j are summation indices over 
elements x, y, z and with 
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- 1 , 


and 


+ l , 
+ 1 , 

- 1 , 


a - + 1 


The terms J Aj are defined as 


and 


J«X - J xx + 1 l COS 2 P 2 • 

jyy = lyy + I £ S i fl 2 R J CO S ^ R J , 

J zz = I ZI + l£ sin 2 p 2 sin 2 Pj , 

Jxy = J xy " l l Sinp 2 COSp 2 COSPj 

Jxz = J xz - l£ sinp 2 cosp 2 sinp, , 

J yz = I yz + If Sin 2 P 2 sinpj COSPj , 


and the elements are 


M + M sin + M cos p 0 ~ ( H H - H H ) - P cos p 7 , 

x b x d 2 z d 1 \ y b z b z b y b / 


S “ M + M x cosp 2 cospj + M sd sin P j - M sinp 2 cosp 1 - (ft H -ft H \ + Psinp 2 cosp x , 

* b d d ' b b b b / 


and 

s 


M z + M x cos p 2 sin p l - M SD cos p x ~ M z ^ si np 2 sin p x - H y ^ H x ^ + P sinp 2 sin^ , 


where P = I £ n Xd ( n y<1 "^2 ] • Since the moment of inertia of the damper rod about its own axis is 
considerably less than its moment of inertia normal to the axis, the I r terms have been dropped. 

It might also be noted that the S b vector terms are still expressed as functions of moments and 
angular velocities referred to the rod coordinate system. This is merely a more convenient way 
of writing out the terms, as representation purely in terms of body coordinates makes the expres- 
sions excessively long. 

The gravitational moment components referred to the frame coordinates and to the rod coordi- 
nates are (from Appendix C): 
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[(I yy -U a 2 a 3 + Ix y a l a 3-I,z a l a 2 +I y z ( a 3 2 - a 2 2 )] ■ 


3 /£ 

»- 3 


[(^“O a l a 3 _I x y a 2 a 3 + I xz ( a l 2 - a 3 2 ) + I y z a l a 2 ] - 


= - 7? [(^x-^y) a l a 2 + I x y ( a 2 2 - a l 2 ) + 1 


xz a 2 a 3 _I yz a 


1 -,] * 


M = 0 , 

X d 


3 M 


My d f 3 l ( b l b 3 ’ 


and 


where 


and 


3/x 

pr bj b 2 , 


a : - cos a cos <fc> , 


a 2 cos a sin^sini// - sin a cos 1 // 


a 3 cos a sin 0 cos 0 + sin a sin 


ij sin p 2 + a 2 cos p 2 cos y o 1 + a 3 cos p 2 sin p^ 


b 2 a 2 sinyOj + a 3 cos p x , 


b 3 - a x cos p 2 ~ a 2 sin/) 2 cospj - a 3 sinp 2 sinpj 


Inspection of Equation 12 reveals that the only unspecified terms remaining are values of p 2 
and p 2 • These values may be determined by considering Equations 6b and 10b, both of which 
have not yet been used: 


1 9 ' J £ n z, M y H + Ms 




n = - n sin^! + n cosp 1 - p 2 

v . V. 1 A 


(6b) 

(10b) 


Elimination of n between the two equations results in 
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( 13 ) 


M y d + M SD 

P.2 = ■ n y b sin ^i + \ cos ^i " n * d n x d ■ r t 

Numerical integration of Equation 13 will produce the term p 2 ; a successive integration then 
yields the value of p 2 , one of the seven system parameters. 

Integration of Equation 12 similarly generates the terms | n x > n v > 1> from which it is 

desired to develop the attitude deviation angles a, </>, and 0 . However, the vector 0 b is referred 
to the inertial frame and may be seen as compounded of both motion of the body relative to the 
attitude coordinate system and motion directly attributable to the motion of the attitude system 
itself. This latter motion must be subtracted from the inertial rotation rate to represent the re- 
lative motion directly related to the attitude deviation angles. The procedure for this is as follows: 

Figure 4b reveals that the angular velocity vectors A, A, and 4> c of the attitude reference 
coordinates with respect to the inertial frame have components in the attitude reference frame of 

£ = A + A sin 4> , 

x o c 

^ y q A + A cos <p c sin A ’ (14) 

- </> sin A + Acos<£ cos A . 

2 o c c 

Since the azimuth angle A is defined as A = arc tan cos 0 c ) , it follows that 

. - (a cos <£ c ) ~ <P C (a cos <p c ~'\<p c sin 4> c ) 

4> 2 + A 2 cos 2 <fi 

c r c 


or 


V A r0 c - v 0 (r A cos <p c - r \<p c sin <£ c ) 


v A 2 + V 


(15) 


In Equations 2 the spherical coordinate functions are expressed in terms of the external forces 
F r ’ f a> F d. • Using these latter terms, Equation 15 becomes 


V A F * - V * F A 


A = 


'(v + v) 


- A sin <fc 


while Equation 14 now takes the form 


(16) 


v A F 0 -v* f a 

c c 

m ( v+ v) 
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(17) 


o f 


and 



cos A + 


sin A 


)• 


The term £ y is equal to 0, since all rotation occurs about the axis from the earth through the 
satellite and the axis normal to the orbit plane. 


With the integrated Equation 12 providing values for fi x , n , and Q z and utilizing Fig- 
ure 5, it follows that 

n = a sin 4> + </> + £ , 

b x b 


0 - a cos 0 sin 0 - 0 cos 0 + £ , 


(18) 


and 


0 - a cos 0 cos 0 + 0sin0 + g 

z b z b 


However, £ b is related to <f 0 the transformation matrix 




(19) 


which yields the equations 

- = ^[( n » b -^ b ) si ^ + ( n *b" f O c ° S0 ] ’ 


and 


* = ' (\ -ty b ) COS * + (\ -^ b ) sin ^ • 

0 = H - £ -a sin 0 , 

x b X b 


( 20 ) 


where all terms on the right hand side will be initially specified. Consequently, numerical inte- 
gration of Equation 20 will produce the attitude deviation angles needed to specify the vehicle 
orientation. 

The damper rod is limited in its plane of rotation by stops placed symmetrically with respect 
to the Y d -axis. When the rod collides with these stops, momentum and energy exchange between 
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the main body and the rod will result. Assuming both total energy and total momentum to be con- 
served, the relative velocity between the main body and the rod after collision is equal in magni- 
tude but opposite in sign to the velocity before collision. This can be shown as follows. Let 


L = I + L COS 2 Pry , 

XX XX t ' 2 

L yy = Iyy + If ( S 1 ^ ^ P t + S itt Vj COS 2 t ) . 

L zz = I.. + If (cos 2 Pj + sin 2 / 3 2 sin 2 Pj) , 

L xy = Iyy " h sinp 2 COsp 2 COSPj , 

L xz = I„ - If sin p 2 cos p 2 sinpj , 


and 


lyz “ If COS 2 p 2 sinpj COS Pj 


Then the total kinetic energy of both bodies is given by 


4 l n 2 + 4l n 2 +4l n 2 +l n n + l n n + l n n 

J XX x k Z yy y b Z zz * b xy x b y b x: x b * b yi y b z h 


+ P 2 If (n yb sinpj - n 2b cospj ) + \ 'p% If 

and the total angular momentum of both bodies is given by 

h\ /° 


/ M 

1 

i 

II 


L L L 

XX x y x z 




L L L 

x y y y y z 


L L L 

y x z y z z 


n 


yb 




Ml 


S in yOj 


y-cos pj 


( 21 ) 


( 22 ) 


Solving Equations 22 for 0,0 , and 0 in terms of p 7 and substituting into Equation 21 

x b y b Z b z 

gives 


p 2 ~ ± y const ant . (23) 

Obviously, one of the values of p 2 given by Equation 23 must be the value when no collision takes 
place. It therefore follows when a collision takes place that 

p 2 (after collision) = -p 2 (before collision) . (24) 
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The values of Ci % , n , and Q z after the collision can then be found by substituting the new value 
of p 2 into Equations (22) and solving simultaneously. 

METHOD OF SOLUTION 

At this point, analysis of the system motion is completed, with the satellite position in space 
determined by integrating Equation (3), the body orientation determined by integrating Equation 
(20), and the angle made by the damper rod determined by integrating Equation (13). The motion 
of the various elements will be continuous and differentiable, except when the damper rod collides 
with one of the stops. This nondifferentiability may be circumvented by utilizing Equations (21) 
and (22) whenever a collision occurs. Appendix D provides a means of finding numerical solutions 
to these equations. 

In the present form the derived equations of motion are functions of general external forces, 
although the gravitational moments derived in Appendix C are considered here to be the only ones 
acting on the system. The development of a general expression allows the inclusion at any time of 
other forces, such as, for example, the thrust from station-keeping thrusters which are studied in 
Part II of this paper. 
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II: A STUDY OF ONE METHOD OF ORBIT 
ADJUST AND STATION KEEPING 


SATELLITE ORBIT AND PARAMETERS 

A medium altitude (approximately 6000 nautical miles) satellite orbit was chosen for the in- 
cluded study. The satellite parameters shown in Table 1 are typical for this orbit.* In Part I of this 
paper a model for the restoring torque, m sd , 
between the main body and the damper rod was 
not given. For this study this torque is as- 
sumed to be: 

m sd = C, P 2 + c 2 p 2 (25) 

C, is a linear viscous damping coefficient and 
C 2 is a linear spring coefficient. It is noted 
that Equation 25 is included in the digital pro- 
gram of Appendix D. 

STEADY STATE DYNAMICS OF THE SATELLITE 


The attitude in which the satellite traverses a circular orbit is shown in Figure 7. This 
motion can best be described as crabbing. The plane of the X, the X b - Y b plane, makes an angle 
of \p ss = -20.76 degrees with the plane of the orbit. The transformation from the body reference 
system to the prime system, included in Figure 7, will be required later. The primed system is 
rigidly attached to the satellite body. 

Attitude librations are induced by eccentricity in the satellite orbit. The amplitudes of these 
librations as a function of eccentricity are shown in Figure 8. The in-plane and cross-plane 
librations correspond to the angles a and $ respectively. The orbit period for each of the data 
points is 6.4 hours. It is noted that for eccentricities of less than 0.05, the increase of in-plane 
amplitude is nearly linear while the cross-plane libration, though initially smaller, is increasing 
at a faster rate. 


Table 1 

Orbit and Satellite Parameters. 


Orbital Period 
Satellite Weight 
I 


if 
p i 

c. 


6.4 hours 
600 pounds 

394.5 slug-ft 2 
3683.0 slug-ft 2 

3288.5 slug-ft 2 
263.1 slug-ft 2 
1.0929 radians 
±1.0 radians 
1.5455 slug-ft 2 /sec 
4.538 slug-ft 2 /sec 2 


*These data were obtained in a private communication with Mr. H. W. Price of the Stabilization and Control Branch at the Goddard Space 
Flight Center in Greenbelt, Maryland. 
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Figure 7— The satellite attitude in a circular orbit 
transformation from the body reference frame to a 
body frame in the orbit plane. 

SATELLITE RESPONSE TO AN IMPULSIVE TORQUE 


Shown in Figure 9 is the response of the satellite to an impulsive torque of 1.3 lb-ft-sec. The 
impulsive torque was applied in the pitch plane and the satellite orbit was circular. It is seen that 
the resulting libration amplitude in pitch is reduced to approximately 10% of its initial value in 
25 hours time. Note that the yaw angle has turned through 180 and is beginning to stabilize near 

i/< ss = - 20.76 - 180.0 = -200.76° . 

This second steady- state yaw angle is possible because of the symmetry of the satellite with re- 
spect to its x b - axis. 


19 



THRUST PROGRAM 


A set of constraints had to be defined before undertaking an orbit adjust and station keeping 
study. In order to take advantage of having one side of the satellite always facing the earth, an 

inversion or turning of the earth pointing axis, 
X b , through 180 degrees must be avoided. The 
determination of a thrust program (thruster 
on-time per pulse and time between pulses) 
required to change the orbit eccentricity by 
0.01 seemed to be a desirable objective. 
Parameters to be optimized should be the 
minimization of the total time required to ad- 
just the eccentricity by the desired amount 
and minimization of the total fuel expenditure. 
Additionally, it would be desirable to obtain a 
thrust program which could be preprogrammed 
and thus eliminate the necessity of ground 
monitoring for successful operation. 

Two one-pound thrusters placed diametri- 
cally opposite each other in direction and with 
respect to the center of mass of the satellite 
were selected for this study.* Each of these 
thrusters was assumed to be 1.5 feet from the 
center of mass. The possibility of either 
tangential or radial thrusting being optimum 
for the above conditions was anticipated, and 
individual studies of each were to be made. For the case of tangential thrusting, one of the 
thrusters was placed on the positive Y b * axis and the other on the negative Y b * axis (see Figure 7). 
For the case of radial thrusting, one of the thrusters was placed on the positive x b axis and the 
other on the negative x b axis. 

The orientation of the thrust vectors, including thrust misalignment, for tangential and radial 
thrusting is shown in Figures 10(a) and 10(b), respectively, with respect to the prime body refer- 
ence frame of Figure 7. The angles and s 2 are defined to give positive components of thrust 
along the positive x b and z b * axes when the angles themselves are positive. For this study a 
total misalignment angle of one degree is assumed. In other words, the magnitudes of S 2 and S 2 
will be such as to combine into a single solid angle of one degree between the thrust vector and the 
reference axis. Combining the definitions of S 2 and s 2 with the transformation of Figure 7 gives 
the forces and moments for tangential thrusting to be: 

F x = T sin Sj cos S 2 , 

b 

•The results of this study can be extrapolated to thrust levels other than one pound if the duration of each pulse is small and the time 
between pulses is the same magnitude when compared to the orbital period; i.e., two pounds of thrust applied for twenty seconds every 
3.2 hours produces approximately the same results as does one pound of thrust applied for forty seconds every 3-2 hours. 



Figure 9-Response of the satellite to an impulsive torque 
of 1 .3 Ib-ft-sec . 
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Figure 10a— Orientation of tangential thruster. 
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and for radial thrusting to be: 


and 



Figure 10b 

T (-cos S, cos 8 2 sin i p ss + sin 
-TP sin S 2 , 

TE sin 8, cos 8- sin d> , 

TE sin 8. cos 8 0 cos </; , 

1 2 r s s 

T cos 8j cos § 2 , 

T (-sin 8j cos 8 2 cos 0 ss + sin 8 2 
T(sin Sj cos S 2 sin i/^ s + sin 8 2 
0 , 

TE (sin S 2 cos i// ss + sin 8j cos 8 2 
TE (-sin 8 2 sin v/> ss + sin Sj cos 8 


—Orientation of radial thruster. 

COS <P ss ) , 


sin «/<„,) , 

cos </< ss ) , 

sin '/'ss) • 

2 COS ^ss) • 



It is noted that the above Equations 26 and 27 are derived for thrusters Located on the negative y’< 
and x b axes, respectively. For sake of simplicity these same equations are assumed for the 
thrusters on the positive Y b > and x b axes, and a negative value, T = -1 pound, is assigned to 
thrust magnitude. Although not shown with the basic digital program in Appendix D, the required 
modification to the program is quite simple. The forces should be evaluated in their body refer- 
ence frame and then transformed to the geocentric coordinate frame. Program locations A(100) 
through A(102) have been reserved for this purpose. The moments should be evaluated in their 
body reference frame and added directly to the program locations A(103) through A(105). 

First order perturbation theory provides a method of determining the most opportune times 
in the orbit for applying thrust pulses to adjust the orbit. To first order the change in eccentricity 
due to tangential and radial perturbation velocities, which can be attributed in this case to the 
thrust pulses, are respectively* 


2 Vi - e 2 cos E y t 
1 - e cos E V » 

(28) 

y 

- y 1 - e 2 sin E . 

E is the eccentric anomaly, e is the eccentricity, V is the satellite velocity, and y t and y r are, 
respectively, tangential and radial perturbation velocities. For small eccentricities, it is seen 
that tangential pulses should be applied at or near apogee or perigee, and that radial pulses should 
be applied at approximately ±90 degrees from apogee or perigee. Thus for the satellite under con- 
sideration, the time between pulses should be one-half the orbital period which is 3.2 hours. It is 
seen also that, all other things being equal, the rate of change of eccentricity for tangential thrust- 
ing is approximately twice that for radial thrusting. 

The one thing remaining in the definition of the problem is to determine the logic required to 
know which of the two thrusters to fire. Since no constraint has been placed on the yaw angle, 
variation between 0 and 360 degrees must be considered for the case of tangential thrusting. 

Table 2 shows the logic required for tangential thrusting to increase eccentricity. The logic re- 
quired for radial thrusting is relatively simple, as the problem was constrained to avoid satellite 
inversion. If the true anomaly of the satellite is between 0 and 180 degrees, the thruster on the 
negative x b axis should be fired, and if the true anomaly is between 180 and 360 degrees, the 
thruster on the positive x b axis should be fired. 


and 


de 

dt 


de 

dt 


*A geometrical discussion of the effects of tangential and radial perturbations is given in Chapter IX of F. R. Moulton’s ** Celestial 
Mechanics *’ , published by The MacMillan Company in 1914. Following similar reasoning, Professor G. D. Boehlor of the Catholic 

University of America has derived Equations 28 in his unpublished lecture notes for Celestial Mechanics. 
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Table 2 


Logic Required for Tangential Thrusting to Increase Eccentricity. \p is 
the Satellite Yaw Angle and 6 is the True Anomaly of the Satellite. The 
Column Under "Thruster” Indicates Whether the Thruster on the Posi- 
tive Y b * Axis or the One on the Negative Y b » Axis Should be Fired. 


\p (degrees) 

d (degrees) 

Thruster 

-110.76 < 0 < + 69.24 

-90 < e <+ 90 

+ 

-110.76 < \p < + 69.24 

+90 < e < + 270 

- 

+ 69.24 <.//<+ 249.24 

-90 < 8 < + 90 

- 

+ 69.24 < xj, < + 249.24 

+90 < e < + 270 

+ 


DISCUSSION OF RESULTS 


Shown in Figure 11 are results of the orbit 
of the graph shows the eccentricity change as a 
time per pulse as parameter. The eccentricity 
a number of computer runs made for the thrust 
alignment angle varying around the reference 
axis (see Figures 10(a) and 10(b)). The points 
are plotted at multiples of 3.2 hours, which is 
one-half the initial orbit period. The circled 
data points indicate tangential thrusting and 
those points enclosed in diamonds indicate 
radial thrusting. The protruding arrows show 
the points which represent a series of runs 
resulting in at least one satellite inversion. 
Examination of the circled points without ar- 
rows reveals that a 0.01 eccentricity change 
can be effected in approximately 40 days with 
either 30-second pulses every 16.0 hours, or 
40-secondpulses every 22.4 hours. The upper 
graph of Figure 1 1 shows that the latter should 
be chosen as less total impulse and, there- 
fore, less fuel is required. Total impulse is 
calculated by simply multiplying the thrust 
level by the total on-time. Also shown on the 
lower graph of Figure 1 1 are two points for 
radial pulses every 22.4 hours. As can be 
seen, 40-second pulses would produce an ec- 
centricity change of 0.0057 in 40 days. This 
latter result was suggested by Equations 28, 


adjust and station keeping study. The lower half 
function of time between pulses, with thruster on- 
changes are average values based on results from 
program indicated, but with the one degree mis- 
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Figure 1 1 -Eccentricity change and total impulse applied 
versus time between pulses. Thruster on-time per pulse 
is shown in parenthesis. The total correction time is 40 
days. 
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after a 0.01 eccentricity change was effected • 
with tangential pulses. 

It is significant that both of the tangential 
data points giving almost a 0.01 eccentricity 
change in 40 days were pulsed at odd multiples 
of one-half the orbital period of 3.2 hours. The 
first multiple was 5 giving 16.0 hours between 
pulses, and the second was 7 giving 22.4 hours 
between pulses. This gives pulses alternating 
between apogee and perigee, which apparently 
has a self-correcting tendency that pulsing al- 
ways at either apogee or perigee does not have, 
as indicated by the lower data points for 18.2 
and 25.6 hours between pulses (see Figure 11). 

Shown in Figure 12 are the increase in 
eccentricity and the drift of the satellite from 
apogee and perigee as a function of time for the 
40-second pulses every 22.4 hours. Here the 
significance is the fact that although the drift 
has increased to approximately 20 degrees, the 
eccentricity is still increasing linearly. This 
means that at least 40 days of preprogrammed 
thrusting can be planned with no penalty being 
paid for drift of apogee and perigee. 

One desirable objective for the tangential 
thrusting program was not achieved. Namely, this investigator did not succeed in finding a method 
of predicting the yaw angle of the spacecraft at the time of the next pulse. (See Table 2 for the 
yaw angle requirements.) No particular pattern seemed to exist in the data. This means that 
ground monitoring or some other method of sensing the yaw angle should be provided for tangential 
thrusting. 

Shown in Table 3 are other pertinent data for the series of computer runs representing 40- 
second pulses every 22.4 hours. It is seen that the number of collisions between the damper rod 
and its stops is quite high for some misalignment angles when pulsing tangentially, while there 
are practically no collisions for radial pulsing. The maximum attitude excursions and the maxi- 
mum attitude angles at the time of impulse are also given. These values seem reasonable. 

CONCLUSIONS 

The results of this study have shown that for the satellite configuration chosen, a thrust pro- 
gram consisting of 40-second pulses every 22.4 hours minimizes the time required to make a 0.01 
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Table 3 


A Summary of the Results of Orbit Correction for Tangential and Radial Thrusting Using 40-second Pulses 
Every 22.4 Hours. The Maximum Values are Maximums for the Entire 40 Day Correction Period. 


S ! 

(Degrees) 

S 2 

(Degrees) 

Maximum a 
(Degrees) 

Maximum <f> 
(Degrees) 

Average 
Number of 
Damper Rod 
Collisions 
Per Pulse 

— 
Maximum 
Collision 
Rate of 
Damper Rod 
(Radians/ 
Second) 

Maximum a 
At Time 
Of Impulse 
(Degrees) 

Maximum <f> 
At Time 
Of Impulse 
(Degrees) 

Tangential Pulses 

0 

1.0 

17.31 

30.47 

1.97 

0.729 x 10 " 3 

5.35 

8.88 

-0.4 

0.9 

31.32 

43.48 

1.62 

1.018 x 10 ' 3 

14.61 

18.09 

-0.6 

0.8 

40.66 

45.11 

1.60 

0.821 x 10 ~ 3 

27.95 

12.07 

-0.7 

0.7 

40.33 

42.30 

0.85 

0.941 x 10 ~ 3 

11.67 

13.80 

-0.8 

0.6 

32.75 

23.41 

0.05 

0.496 x 10 ' 3 

3.11 

2.41 

-0.9 

0.4 

66.90 

46.38 

0.45 

0.523 x 10 -3 

18.73 

13.68 

-1.0 

0 

53.74 

33.05 

0 

— 

10.17 

12.06 

Radial Pulses 



' 




0 

1.0 

36.44 

37.61 

0 



5.89 

12.29 

-0.4 

0.9 

56.41 

51.17 

0.03 

0.181 x 10 ~ 3 

23.24 

17.81 

-0.6 

0.8 

42.16 

39.28 

0 

— 

9.16 

13.95 

-0.7 

0.7 

40.09 

41.71 

0 

— 

6.57 

15.14 

-0.8 

0.6 

46.23 

42.87 

0 

— 

8.05 

9.47 

-0.9 

0.4 

44.18 

26.29 

0 

— 

8.72 

7.95 

-1.0 

0 

55.36 

36.12 

0.02 

0.099 x 10 ~ 3 

9.83 

6.93 


eccentricity change in the orbit and minimizes the fuel required. A one pound thrust level and a 
one degree thruster misalignment angle were assumed for the study. If tangential thrusting is 
used, the 0.01 eccentricity change can be effected in 40 days. However, a means of sensing the 
yaw angle of the satellite must be provided in order to establish the direction in which successive 
pulses should be made. Additionally, tangential thrusting results in a significant number of col- 
lisions between the damper rod and its mechanical stops. No attempt was made to evaluate these 
collisions and rates of collision in terms of material strengths. 

Should sensing the yaw angle and the damper rod collisions prove to be unsatisfactory, radial 
thrusting can be used. Both of the problems, sensing the yaw angle and the rod collisions, are 
avoided, but at the expense of almost doubling the time required to adjust the orbit. In 40 days, 
radial thrusting produces an eccentricity change of only 0.0057. 


(Manuscript received February 9, 1966) 
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Appendix A 

List of Symbols 


A 

A 

Ci 

c 2 

F 

E 

e 

H 

h t 

I 

U. j, k} 

J 

K 

L 

M 

m sd 

m 

R 

{ r - A - ^c} 

T 

w 

V 

{X, Y, Z} 
{a, 4>, i p) 


acceleration vector of satellite center of mass 
azimuth angle 

linear viscous damping coefficient 

linear spring coefficient 

external forces acting on the satellite 

eccentric anomaly of the satellite 

eccentricity of the satellite orbit 

total angular momentum of the satellite main body 

total angular momentum of both the main body and the damper rod 

moment of inertia tensor of satellite main body 

unit vectors along coordinate axes X , Y , z 

total moment of inertia tensor including damper rod contribution 
inverse of the total moment of inertia tensor 

the inertia tensor for momentum transfer between the two satellite sections at the 
time of collision 

gravitational moment exerted on the satellite by the earth 

restoring torque exerted on the rod by the spring-damper mechanism 

mass 

reaction torque between the two satellite sections 
inertial spherical coordinates of the satellite 
thrust magnitude of orbit adjust thrusters 

matrix defining vector transformation from X -base coordinate system to ! -base 
system 

inertial velocity vector of the satellite center of mass 
orthogonal right-handed coordinate axes 

pitch, roll, yaw angles of satellite measured from attitude reference system 
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steady state yaw angle 

determinant of the total moment of inertia tensor 

thrust misalignment angles 

tangential and radial perturbation velocities 

total external moments acting on the satellite 

earth gravitational constant 

inertial angular velocity vector of the satellite 

dihedral angle between rod rotation plane and satellite X plane 

angular deviation of damper rod from equilibrium position 

maximum allowable deviation angles made by rod before hitting stops 

true anomaly of the satellite 

Subscripts 

refers to main body coordinate system 

refers to geocentric coordinate system 

refers to damper rod coordinate system 

refers to attitude reference or orbital coordinate system 

refers to the length of the damper rod 

refers to the radius of the damper rod 

refers to inertial spherical coordinates of the satellite 

refers to a specified coordinate system 


28 


b 

c 

d 

o 

Z 


{r. A, </> c } 
{x, y, z} 


0 

r ss 

A 

( r t- n} 

i 

H 

Pi 

P2 

~Ps 

e 
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Appendix B 

Elements Defining a Satellite Orbit 


The orbit of a satellite around the earth is, according to Kepler's first law, an ellipse with 
the earth located at one focus. The semimajor axis (a) and the eccentricity (e) define the size 
and shape of the ellipse, while the parameter (r) denotes the time of perigee passage of the 
satellite. These three elements permit calculation of the vehicle position in its orbit plane at 
any time. 

An additional set of three variables is also required to orient the orbit plane itself. As shown 
in Figure Bl, the longitude of the ascending node (n) determines the line of intersection of the orbital 
plane and the earth's equatorial plane, the element (I) specifies the angle at which the orbital plane 
is inclined to the equator, and the longitude of perigee (w + 0) determines the direction of perigee 
with respect to the vernal equinox. The six elements a, e, r, fi, I, and (a> + 6) are called orbital 
elements, and all six are required to completely specify the satellite orbit in relation to the earth. 


These elements designate position as a func- 
tion of time, and therefore implicitly contain the 
satellite velocity components. The solution of the 
translational equations of motion (Equation 3) re- 
quires velocity- dependent and position- dependent 
terms as initial conditions, but referred to the 
inertial coordinate system. Conversion from 
orbital elements to inertial elements may be ef- 
fected as follows: 

(1) Position elements r, <p c , A : 

For r, 

_ a(l~ e 2 ) 

r 1 + e cos 6 (Bl) 

(0 = true anomaly; determined 
from e and r); 


for 0 c> - 


sin 0 

r c 


sin I sin (co + 0) 


(B2-a) 



cos 4> c > 0 ; 


Figure Bl —Rotational transformation from the orbit 
or attitude to the inertial reference frame. 


(B2-b) 
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for A, 


sin A = cos cp [cos cos I sin (a) + 0) + sin H cos (a) + 0)] 


(B3-a) 


cos A = cos ^ [cos ft cos (ai + 0) - sin H cos I sin (co + 0)J 


(B3-b) 


(2) Velocity elements r, r 2 4> , r 2 A cos 2 4>, * 


For r, 


(Reference 7, p. IE-22). The two remaining terms arise in the following manner: Noting that from 
Figure 4b: 


A cos 4> c ~ € z cos A , 


(B5-a) 


0„ - sin A , 


(B5-b) 


and from Reference 7, p. ni-23, 


« 2 f. 




( 1 + e cos 0) 2 


[fjL( 1 + e cos 6) 


Solution for the two remaining constants by substitution yields, for r 2 <f> : 


for r 2 A cos 2 cp : 


Vr(l + e cos 6 ) sin A 


r 2 Acos 2 <£ - i J jLtr( 1 + e cos 0) cos c£ c cos A . 


After internal calculations have been completed in inertial coordinates, conversion to orbital 
elements may be done through the following equations: 

For I: 

cos I = cos <b cos A , 

’ c 




For n : 


For (a) + 8): 


sinft - (sin A sin A - cos A sin <£ c cos A )/sin I , 
cos = (cos A sin A + sin A sin cos A)/sin I ; 


sin(ciJ + 0) - sin<£ c /sin I , 

cos(a> + 0) = (cos 0 C sin A)/sin I . 


When the orbit plane coincides with the earth equatorial plane (I = 0), then <t> c = A = o. 
ing the longitude of ascending node (n)to be zero, it follows that 


a + 6 = A . 

Defining the element Q = rv 2 /u and cos 2 r = V A 2 -V ^ 2 j V 2 , the determination of elements 
are: 


a - 2 - Q ’ 


and 


e = yi-Q(2-Q) cos 2 T 

If the orbit is not circular (e / 0), it may be shown that 



If the orbit is circular, define u = 0 . 


(BIO) 

(Bll) 
By defin- 

(B12) 

(a) and (e) 

(B13) 

(B14) 

(B15) 
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Appendix C 

Gravitational Moments Generated by Earth Potential 


Roberson (References 8 and 9) has shown the rotational potential of a satellite vehicle under 
the influence of a y,f r 2 gravity force field to be 


U 


3 mI 0 
Tr 1 ’ 


(Cl) 


where i 0 is the instantaneous moment of inertia of the satellite about the local geocentric vertical. 
When expressed in terms of the body moments of inertia and the attitude deviation angles a, 0, and 
0 , the result is 


Ixx*l 


+ I a~ 
yy 2 


+ I 


21 a 
*y 


i d 2 


2I xz a i a 3 


+ 21 


1 2 3 


(C2) 


where 


a l - cos a cos 0 , 

a 2 - - cos a sin 0 sin 0 - sin a. cos 0 , 

and 

a 3 = - cos a sin 0 cos 0 + sin a sin 0 . 

DeBra (Reference 10) points out that work is related to a change in potential by 

M * da T = - dU , 

where M = total moment, and 

da T = resultant of incremental changes in da, d0, and dp. 

Expressing the elements da, d0, d0 in body coordinates gives 

|"«9U <9U"| 

M b * B • da T - " [da ’ d<p ' ^0 J ’ da r ’ 


(C3) 


(C4) 
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or 


M b • B = 


where M b denotes the total moment, now referred to body coordinates; and the dyadic B is given by 
the matrix, 


B - I cos 0 sin 0 - cos 0 0 

\ cos 0 cos 0 sim p 0 


It follows from Equation C-4 that 


5U sin 0 dU <9U sin 0 sin 

da cos 0 + <90 COS( P + <90 cos 0 


dV cos 0 dl) m <9U sin 0 cos 0 

da cos <56 ” 50 n r + 50" cos 0 


which gives the expressions 


[( J yy - 1 zz) a 2 a 3 + J xy a l a 3 " a l a 2 + l yz { S 


r 3 [(^11 O a l a 3 *xy a 2 a 3 + *xz ( a l 2 3 3 ) + *yz 3 1 


_ ~ 


~Jt [(^zx" X yy) a l a 2 + I xy ( a 2 2_a l 2 ) + I x 


xz a 2 a 3" I yz a l 
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The moment referred to the rod coordinates may be calculated by noting that 




} d 
k d 





The moment terms, when referred to this system, are given by 

M = 0 , 

x d 


I'd 


b, b. 


7T l£ 


where 


*>i 

b 2 

b 3 


a j sin p 2 + a 2 cos p 2 cos p x 
- ajsin/Oj + a 3 cosp 1 , 
aj cos p 2 - a 2 sin p 2 cos p x 


+ a 3 cos p 2 sin p x , 


- a 3 sin p 2 sin p x . 


(C9) 
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Appendix D 


A Gravity-Oriented Satellite Digital Program 


INTRODUCTION 

The motion of a two-body, gravity- oriented satellite has been implicitly determined by a set 
of seven second-order differential equations. Three of these equations locate the satellite in an 
inertial reference system, three specify its attitude with respect to a designated attitude reference 
system, and the seventh describes the motion of the damper rod in relation to the main body. 

A numerical technique which facilitates obtaining useful results is now presented. The pro- 
gram, developed along lines similar to a program by Dennison and Butler*, is written in FORTRAN 
n format and occupies approximately 7000 words of memory on the SDS-920 computer. This in- 
cludes all necessary standard library subroutines, as well as the main program and attendant 
subroutines. 

A listing of the input variables, to determine the initial conditions of the equations, and vehicle 
parameters is included. From these the single-step Runge-Kutta quadrature program generates 
successive values of the orbit location and vehicle attitude angles, along with the rod deviation 
angle. A choice of output options is provided. 

The program as it presently stands includes only gravitational forces and moments generated 
by the earth on the satellite, but it can be readily adapted to accommodate additional forces and 
moments from other sources through the addition of an intermediate subroutine. To facilitate this 
or any other change that might be desired, internal variables are listed and defined, and the pro- 
gram listing (see Table DIO) is interlaced with comment statements. 

EQUATION LISTING 

The previously derived equations defining the satellite motion are described briefly in this 
section. Reference should be made to the derivation for secondary equations, and to the accompa- 
nying list of symbols for definitions. 


♦Dennison, A. J., and Butler, J. F., "Missile and Satellite Systems Program for the IBM-7090**, Technical Report 61-SD-170 (The Gen- 
eral Electric Co., Missile and Space Division, King of Prussia, Pa.), prepared on contract to NASA, unpublished. 
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Translational Motion 


The variables (r, A, 4> c ) are determined by doubly integrating Equation(s) 3, which are re- 
peated here as: 

F = (r A 2 cos 2 4 > c + r<£ c 2 ) + -jj- ' 



cos 2 <p c ) 



r cos 


*.) * 


_d_ 

dt 


(r’ij 



r A 2 sin cos <p 


(Dl) 


Rotational Motion 

The inertial body rate, fi, is determined by integrating Equation (12), which is repeated here as: 


K\ 


( K K K \ 

' xx xy x: ] 

M 


n 

^b 

i 

- A 

K K K 

xy yy yz 

i 

a 

Vb 

(D2) 

\ h J 


\K K K / 

\ x z yz z z / 

V-j 



The attitude deviation angles are then determined by integrating Equation(s) (20), which are re- 
peated here as: 


“ = ^b*[( n y b “O sin ^ + ( n * y -^ y ) cos ^] ’ 
* = “ K -*>*) “•t+K-Q sinxp ’ 

\b = Q - A - a sin <£ . 

x b x b 

Rod Motion 


(D3) 


The rod motion, defined by the variable p 2J is found by doubly integrating Equation (13), which 
is repeated here as: 


f>2 


- n 


*b 


sin p 1 


+ COS Pi 

z b 


- n n 



(D4) 
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In Equation (D4), the spring restoring torque M SD is defined by the equation: 

Mgp - Cj f>2 + C 2 P2 • 

Since p 2 is dimensionally a pure numeric, it 
follows that c, has the dimensions 


[m] M [r 1 ] , 

and C 2 has the dimensions 

[ m ] U 2 ] [t- 2 i . 

COMPUTER PROGRAM 
Flow Chart 

The main program and associated subrou- 
tines are related, as shown in Figure Dl. 

Program Description 

The Main Program [ MAIN ] performs three 
major functions prior to calling the integration 
subroutine: (a) initial read- in and/or computa- 
tion of all constant vehicle parameters used in 
the program, values stored in the (A) matrix 
from A(50) through A(97), (b) read-in of initial 
conditions for rotational equations of motion 
and for equation of damped- rod motion, and (c) 
read-in of orbital elements and conversion to 
initial conditions in the inertial coordinate sys- 
tem for the translational equations of motion. 



Figure Dl— Main program and associated subroutines for 
numerical determination of the equations of motion of a 
gravity-gradient satellite. Standard FORTRAN II library 
subroutines used are arctangent, cosine, sine, and square 
root. 


The vehicle parameter input format is listed in Table Dl, and is entered into the computer 
after WRITE OUTPUT TAPE 3, 2 of the main program. 

The input format for initial conditions and orbital elements is listed in Table D2, and is 
entered into the computer after A(97) = A(62)**2 of the main program. Following the evaluation 
of these inputs, control is transferred to the integration subroutine RK(N) . 
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Table D1 

Vehicle Parameter Input Format. 


ENTRY TAPE 

ENTRY CODE 

PARAMETER 

UNITS OF MEASURE 

READ INPUT 

A(50) 

Vehicle mass 

Slugs 

TAPE 2,1, 

A(51) 

Rod mass 

Slugs 

READ INPUT 

A(53) 

Vehicle inertia moment (x) 

SIug-ft. 2 

TAPE 2,1, 

A(54) 

Vehicle inertia moment (y) 

Slug-ft. 2 


A(55) 

Vehicle inertia moment (z) 

Slug-ft. 2 

READ INPUT 

A(59) 

Vehicle inertia product (xy) 

Slug-ft. 2 

TAPE 2,1, 

A(60) 

Vehicle inertia product (xz) 

Slug-ft. 2 


A(61) 

Vehicle inertia product (yz) 

Slug-ft. 2 

READ INPUT 

A(62) 

Rod inertia moment 

Slug-ft. 2 

TAPE 2,1, 

A(63) 

Damper plane inclination 

Degrees 

READ INPUT 

A(66) 

Cl 

Slug-ft. 2 /sec. 

TAPE 2,1, 

A(67) 

c 2 

Slug-ft. 2 /sec. 2 


Table D2 

Input Format for Initial Conditions and Orbital Elements. 


ENTRY TAPE 

ENTRY CODE 

PARAMETER 

UNITS OF MEASURE 

READ INPUT 

AA 

Semimajor axis 

Feet 

TAPE 2,1, 

ECC 

Eccentric ity 

Numeric 


XI 

Orbit Plane Inclination 

Degrees 

READ INPUT 

THE 

True Anomaly 

Degrees 

TAPE 2,1, 

AP 

Argument of Perigee 

Degrees 


AN 

Longitude of Ascending Node 

Degrees 

READ INPUT 

Y(7) 

Inertial body rate (x) 

Rad. /sec. 

TAPE 2,1, 

Y(8) 

Inertial body rate (y) 

Rad. /sec. 


Y(9) 

Inertial body rate (z) 

Rad. /sec. 

READ INPUT 

Y(10) 

Pitch angle 

Degrees 

TAPE 2,1, 

Y(ll) 

Roll angle 

Degrees 


Y(12) 

Yaw angle 

Degrees 

READ INPUT 

Y(13) 

Rod angular velocity 

Rad. /sec. 

TAPE 2,1, 

Y(14) 

Damper rod angle 

Degrees 

READ INPUT 

Y(15) 

Machine integration time 

Seconds 

TAPE 2,1, 

Y(16) 

Integration interval 

Seconds 


Y(17) 

Printout interval 

Seconds 


Y(18) 

Termination time 

Seconds 


The Runge-Kutta Subroutine [RK(N)] performs the fourth-order integration procedure. Upon 
completion of this subroutine function, control is transferred back to the main program, resulting 
in the termination of the computer program. 

Entry into the Derivative Subroutine [. DERIV ] is from subroutine RK(N) for the purpose of 
computing the values of the differential equations at each successive point for each differential 
equation. Control is then returned to subroutine RK(N). 
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. Entry into the Momentum Transfer Subroutine [. AUTO (FLAG, STP)] is from subroutine 
RK(N), for the purpose of computing momentum transfer when the damper rod collides with the 
stops mounted on the main vehicle. Time of collision and relative angular velocity between the 
sections are written on output tape; control is then returned to subroutine RK(N). 

The Data Output Subroutine [OUT(K)] is called by subroutine RK(N) in order to convert inte- 
grated values of the sought-after variables to the proper output format. The translational coordi- 
nates are converted to orbital element representation, while all angular measures are converted 
from radians to degrees. Depending on the value of K (either 0 or 1), one of two possible output 
modes, described in Table D3, is activated. 


Table D3 


Data Output Subroutine. 


OUTPUT TAPE 

CODE 

PARAMETER 

UNITS 

OUTPUT MODE (0) 

WRITE OUTPUT 

THR 

Time 

Hours 

TAPE 3,32, 

PITCH 

Pitch angle 

Degrees 

(LINE 30) 

ROLL 

Roll angle 

Degrees 


YAW 

Yaw angle 

Degrees 


RH02 

Damper rod angle 

Degrees 

OUTPUT MODE (1) 

WRITE OUTPUT 

THR 

Time 

Hours 

TAPE 3,20, 

Y(15) 

Time 

Seconds 

WRITE OUTPUT 

PITCH 

Pitch angle 

Degrees 

TAPE 3,21, 

Y<7) 

Inertial body rate (x) 

Rad. /sec. 


AA 

Semimajor axis 

Feet 


THE 

True anomaly 

Degrees 

WRITE OUTPUT 

ROLL 

Roll angle 

Degrees 

TAPE 3,21, 

Y(8) 

Inertial body rate (y) 

Rad. /sec. 


ECC 

Eccentricity 

Numeric 


AP 

Argument of perigee 

Degrees 

WRITE OUTPUT 

YAW 

Yaw angle 

Degrees 

TAPE 3,21, 

Y(9) 

Inertial body rate (Z) 

Rad. /sec. 


XI 

Orbit inclination 

Degrees 


AN 

Longitude of ascending node 

Degrees 

WRITE OUTPUT 

Y(4) 

Radius 

Feet 

TAPE 3,20 

V 

Velocity 

Ft. /sec. 


RH02 

Damper rod angle 

Degrees 


Y(13) 

Rod angular velocity 

Rad. /sec. 


MATRIX ELEMENTS 

Most of the physical quantities computed in the program are assigned to one- dimensional ar- 
rays "A” and "Y". This eliminates the necessity of inventing meaningful alphabetic characters to 
name each of these quantities individually. A less obvious advantage to persons familiar with 
FORTRAN is the ease with which information can be transferred between subroutines by simply 
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pu tting these arrays in "COMMON". It is noted that several of the array elements are intentionally 
undefined in the present program, to facilitate addition of features such as special thrust routines, 
higher order gravitational calculations, and a more exact mathematical model of the spring-damper 
mechanism. 

The first fourteen elements of the Y-array have been assigned to the basic variables for the 
seven-degree-of-freedom equations of motion. The time derivatives of these variables have been 
assigned to the corresponding A-array elements indicated in Table D4. 


Table D4 

Elements of Translational and Rotational Motion. 


ELEMENT 

PARAMETER 

ELEMENT 

PARAMETER 

A( 1) 

r 

Y( 1) 

r 

A( 2) 

3T (r 2 A cos 2 4> c ) 

Y( 2) 

r 2 A cos 2 <p 

c 

A( 3) 

3t ( r2 ^c) 

Y( 3) 

r 2 k 

A( 4) 

r 

Y( 4) 

r 

A( 5) 

A 

Y( 5) 

A 

A( 6) 


Y( 6) 


A( 7) 

n 

x b 

Y( 7) 

n 

x b 

A( 8) 

n 

Y( 8) 

n 

yb 

A( 9) 

n 

z b 

Y( 9) 

n 

z b 

A(10) 

a 

Y(10) 

a 

A(H) 

4> 

Y(ll) 

4> 

A(12) 

+ 

Y(12) 

0 

A( 13) 

P2 

Y(13) 

P2 

A(14) 

P2 

Y(14) 

P2 


Integration Control Parameters 

The integration control parameters are as follows: 

Y(15) = Current machine integration time, 

Y(16) = Integration time interval, 

Y(17) = Printout time interval, 

Y(18) = Termination time. 
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Satellite Inertial Parameters 


The satellite inertial parameters A(50) through A(97) in Table D5 remain constant for any 
given run. The repeated unnecessary entry into the sine-cosine subroutine and other arithmetic 
operations required for A(75) through A(97) is avoided by making these calculations early in the 
program (Main Program) and retaining them for later use. The resulting time saved is 
appreciable. 


Table D5 

Vehicle Inertia Parameters. 


ELEMENT 

PARAMETER 

ELEMENT 

PARAMETER 

A(50) 

m b 

A(75) 

ii - 1 2 

yy zz yz 

A(51) 

m d 

A(76) 

II - I I 

xz yz zz xy 

A(52) 

m b + m d 

A(77) 

^xy ^yz ^yy ^ xz 

A(53) 

!xx 

A(78) 

( I + la ) I - I 

\ XX t / ZZ XZ 

A(54) 

*yy 

A(79) 

!* y Ixz - (Ixx +I f) !yz 

A(55) 

I,. 

A(80) 

(Ixx + If) !yy- Ixy i 

A(56) 

J yy- 

A(81) 

If (l yy s in 2 Pi - 2 I yi sin p, cos p x discos 2 p,) 

A(57) 

Ixx 

A(82) 

If (-I y2 sinp! +I„cosp 1 ) 

A(58) 

^xx *yy 

A(83) 

If (" Ixy sin Pi + I XI COS Pi) sinp! 

A(59) 

I 

xy 

A(84) 

If (lyy Sin /°l ' lyx COS Pl) 

A(60) 

I 

xz 

A(85) 

If (ixy Sin Pi ~ X xx OOSPj) COSPl 

A(61) 

!y. 

A(86) 

2 If Ixx Sin ^l 

A(62) 

If 

A(87) 1 

If (lxx sin2 ^l "Ixx) 

A(63) 

p 1 

A(88) 

If (-Ixx Sin Pl COS Pl + lyx) 

A(64) | 

sin pj 

A(89) 

" !f ( J xy Sin /°1 +I xx COS Pl) 

A(65) 

COS p x 

A(90) 

2 ^ X xy CO *Pl 

A(66) 

c, 

A(91) 

If (lxx Cos2 Pl -lyy) 

A(67) 

C 2 

A(92) 

If (l yy cos 2 Pi + 2 I yz sinpjcospj +I„sin 2 p,) 

A(68) 

I + Ifl 

XX t 

A(93) 

- If (l xy cosp, + I„ Sinp,) 

A(69) 

sin 2 p j 

A(94) 

If r xx cos2 -°l 

A(70) 

2 

cos p x 

A(95) 

If X xx Sitl Pi COS Pi 

A(71) 

sin p x cos p 1 

A(96) 

If I.xx sir,2 £l 



A(97) 

If 2 


Translational Equation Forces 

The translational equation forces are as follows: 
A(100) = F r , 

A(101) = F a , 
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A(102) = . 

The program described herein contains zeroes for these elements. Any type of force (lunar and 
solar perturbations, for example), which the user might desire, should be added to these elements. 
These forces should be evaluated in Subroutine Deriv. 

Rotational Equation Moments and Variables 

The rotational equation moments and related variables are given in Table D6. The moments 
are assigned to elements A(103) through A(109). The other variables are angular momentum, 
A(110)-A(112) and A(116); angular velocities, A(113)-A(115) and A(117)-A(125); damper rod position 
functions A(126) and A(127); and the inertia tensor A(128)-A(134). These are required to evaluate 
the moments. These calculations are made in Subroutine Deriv . 

Momentum Transfer Terms for Collision 

Table D7 relates momentum transfer terms for collisions between the vehicle body and the 
damper rod. This evaluation is made in Subroutine Auto. 


Table D6 

Rotational Equation Moments and Variables. 


ELEMENT 

PARAMETER 

ELEMENT 

PARAMETER 

ELEMENT 

PARAMETER 

A(103) 

M 

x b 

A(115) 

n 

z a 

A(126) 

sin 2 p 2 

A(104) 

M 

y b 

A(116) 

i„ n n 

1 x d Vd 

A(127) 

sin p 2 cos p 2 

A(105) 

M 

2 b 


" 1 i ^x d 2 

A(128) 

K x/ A 

A(106) 

M 

x d 

A(117) 

s 

x b 

A(129) 

K /A 

xy/ 

A(107) 

M 

y d 

A(118) 

s 

y b 

A(130) 


A(108) 

M 

z d 

A(119) 

8 

*b 

A(131) 

K y /A 

A(109) 

m sd 

A(120) 

^ x o 

A(132) 


A(110) 

H 

x b 

A(121) 


A(133) 

k 2 /a 

A(lll) 

H 

y b 

A(122) 


A( 134) 

A 

A(112) 

H 

2 b 

A(123) 

CO 

x b 



A(113) 

n 

x d 

A(124) 

CO 

y b 



A(114) 

n 

y d 

A(125 





44 








Table D7 

Momentum Transfer Terms. 


ELEMENT 

PARAMETER 

ELEMENT 

PARAMETER 

A(135) 

h t 

X 

A(140) 

L 

xz 

A(136) 

H t + (pi) l h sin^j 
y 

A(141) 

L 

y 

A(137) 

H T z - (- 02 ) 1 l t MS Pi 

A(142) 

L 

yz 

A(138) 

L 

X 

A(143) 

L 

z 

A(139) 

L 

xy 

A(144) 

V 


Trigonometric Functions 

Table D8 groups the trigonometric functions required by the program. Elements A(226) 
through A(231) are combinations of the previously defined trigonometric functions which are re- 
quired for the gravitation moments. 

Program Listing and Symbol Table 

Table D9 defines alphabetic symbols which, for the most part, are used internally in the pro- 
gram. The FORTRAN n digital program listing is shown as Table DIO. 


Table D8 

Trigonometric Functions. 


ELEMENT 

PARAMETER 

ELEMENT 

PARAMETER 

ELEMENT 

PARAMETER 

A(200) 

sin A 

A(201) 

cos A 

A(220) 

cos cp sin \p 

A(202) 

sin <p 

r c 

A(203) 

COS0 c 

A(221) 

cos < p cos l p 

A(204) 

sin A 

A(205) 

cos A 

A(222) 

cos p 2 cos p 1 

A(206) 

sin a 

A(207) 

cos a 

A(223) 

cos p 2 sin p x 

A(208) 

sin <p 

A(209) 

cos cp 

A(224) 

- sin p 2 cos p x 

A(210) 

sin 0 

A(211) 

cos i// 

A(225) 

- sin p 2 sin p x 

A(212) 

sin p 2 

A(213) 

cos p 2 

A(226) 

a l a 2 

A(214) 

cos a cos cp 



A(227) 

a l a 3 

A(215) 

- cos a sin 0 sin 1 // - sinacosi/i 



A(228) 

a 2 a 3 

A(216) 

- cos a sin 0 cos 0 + sinasin \p 



A(229) 

a \ 

A(217) 

sin a cos <p 



A(230) 

a i 

A(218) 

- sin a sin 0 sin 0 + cos a cos i// 



A(231) 

*1 

A(219) 

- sin a sine/) cos \p - cos a sin ^ 
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Table D9 

Alphabetic Symbols. 


AA semimajor axis of the satellite orbit 

AN longitude of the ascending node of the satellite orbit 

AP argument of perigee of the satellite orbit 

CAN cosine of the longitude of the ascending node 

CAP cosine of the argument of perigee 

CTHE cosine of the true anomaly 

CXI cosine of the orbit plane inclination 

ECC eccentricity of the satellite orbit 

FLAG a storage location used to transmit the signal of a rod collision with the main body 
P semilatus rectum of the satellite orbit 

PITCH attitude angle of the satellite 

RAD conversion factor for converting radians to degrees 
RH02 angular deviation of the damper rod from equilibrium position 
ROLL attitude angle of the satellite 

R2 radius vector squared 

R3 earth constant divided by radius vector cubed 

SAN sine of the longitude of the ascending node 

SAP sine of the argument of perigee 

STHE sine of the true anomaly 

STP indicator of the status of calculation of momentum transfer during damper rod collision with the 
main body 

SXI sine of the orbit plane inclination 

THE true anomaly of the. satellite in its orbit 

THR time in hours 

V total velocity 

VL longitudinal component of velocity 

VP latitudinal component of velocity 

VT vector sum of longitudinal and latitudinal components of velocity 

XI orbit plane inclination 

YAW attitude angle of the satellite 

ZMU earth's gravitation constant 

ZMU3 minus three times the earth's gravitation constant 
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Table DIO 
Program Listing. 


C MAIN PROORAM 

DIMENSION A( 500) . Y( 100) 

COMMON A » Y » R AD » ZML- » 2MU3 
RAD = 57.295779 
ZMU = 1.40V6449E16 
ZMU3 = -3 • *ZMU 

1 FORMAT ( At 13 « 8 ) 

2 FORMAT (1H1) 

C VtHICLt PARAMETER INPUT (TABLE D1 ) 

WRITE OUTPUT TAPE 3.2 

READ INPUT TAPE 2 , 1 > A < 50 ) . A ( 51 ) 

READ INPUT TAPE 2 ♦ 1 . A ( 5 3 ) » A ( 54 ) , A < 55 ) 

READ INPUT TAPE 2 . 1 * A ( 59 ) * A < 60 ) ♦ A ( 61 ) 

READ INPUT TAPE 2 , 1 * A ( 62 ) * A ( 63 ) 

READ INPUT TAPE 2 . 1 . A< 66 ) > A ( 67 ) 

WRITE OUTPUT TAPE 3. 1 * A ( 50 ) ♦ A ( 51 J 
WRITE OUTPUT TAPE 3 * 1 » A ( 53 ) > A ( 54 ) » A ( 5 5 ) 

WRITE OUTPUT TAPE 3 * 1 » A ( 59 ) * A ( 60 ) , A ( 6 1 ) 

WRITE OUTPUT TAPE 3 * 1 * A ( 62 ) > A < 63 ) 

WRITE OUTPUT TAPE 3 » 1 » A ( 66 ) » A < 67 ) 

C CALCULATION UF CONSTANT VEHICLE INERTIA PARAMETERS (TABLE D5 ) 

A ( 52 ) = A ( 50 ) +A (51) 

A ( 56 ) = A ( 54 ) -A ( 55 ) 

A ( 57 ) = A I 55 ) -A ( 53 ) 

A ( 58 ) = A ( 5 3 ) -A ( 54 ) 

A l 65 ) = A ( 63 ) /RAD 
A ( 64 ) = S INF ( A ( 65 ) ) 

A (65) = COSF ( A ( 6 5 ) ) 

A ( 68 ) = A ( 53 ) +A ( 62 ) 

A ( 69 ) = A(6A)**2 
A ( 70 ) = A ( 65 ) **2 
A ( 7 1 ) = A ( 64 ) *A ( 65 ) 

A ( 7 5 > = A(54)*A(55)-A(61 )**2 
A 1 76 ) = A(60)*A161)-A(55)*A(59) 

A ( 77 ) = A(59)*A<61)— A(54)*A(6Q) 

A ( 78 ) = A(68)*A(55)-A<60)**2 
A ( 79 ) = A(59)*A(60)~A(68 ) *A ( 61 ) 

A l 80 ) = A(68)*A(54)-A(59)**2 

A ( 8 1 ) = A(62)*IA(54)*A(69)~2.*A(61)*A(71)+A(55)*A(70)) 

A ( 82 ) = A(62)*(A(55)*A(65) -A ( 61 ) * A ( 64 ) ) 

A ( 8 5 ) = A ( 62 ) * ( A( 59 )*A(64)— A(60)*A(65 ) ) 

A ( 8 3 ) = -A ( 85 ) * A ( 64 ) 

A ( 84 ) = A(62 ) * 1 A ( 54 ) *A l 64)-A(61)*A(65 ) ) 

A ( 8 5 ) = A ( 85 ) *A ( 65 ) 

A ( 86 ) = 2.*A(62)*A(60)*A<64) 

A ( 87 ) = A(62)*(A(53)*A(69)-A(55) ) 

A ( 88 ) = A ( 62 ) * ( -A ( 53 ) *A ( 7 1 ) +A( 61 ) ) 

A l 89 ) = -A(62)*(A(59)*A(64)+A(60)*A(65) ) 

A ( 90 ) = 2.*A(62)*A(59)*A(65) 

A ( 91 ) = A(62)*(A(53)*A(70)~A(54) ) 

A ( 92 ) = A(62)*(A<54)*A<70)+2.*A<61)*A(71)+A(55)*A(69) ) 

A ( 9 3 ) = — A'( 62 ) * ( A ( 59 ) *A ( 65 ) + A ( 60 ) *A ( 64 ) ) 

A ( 96 ) = A ( 53 ) *A ( 62 ) 

A ( 94 ) = A ( 96 ) *A ( 70 ) 

A ( 95 ) = A ( 96 ) *A (71) 

A l 96 ) = A ( 96 ) *A ( 69 ) 

A ( 97 ) = A ( 62 ) **2 

INITIAL CONDITION ANu uKoiTAL ELEMENT INPUT (TABLES U2 AND D4 ) 
Y ( 1 ) THRU Y ( 14 ) ARE ALL INITIALIZED 
READ INPUT TAPE 2 , 1 » AA , ECC »X I 
READ INPUT TAPE 2 , 1 . THE » AP * AN 
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Table DIO (Continued) 
Program Listing. 


READ INPUT TAPb 2 » 1 * Y { 7 ) 9 Y ( 8 ) 9 V ( 9 ) 

READ INPUT TAPb 2 * 1 9 Y ( 10 ) 9 Y ( 1 1 ) * Y ( 12 ) 

READ INPUT TAPE 2 * 1 » Y ( 1 3 ) ♦ Y ( 14 ) 

READ INPUT TAPb 2 * 1 * Y < 1 5 ) * Y < 16 ) t Y ( 17 ) » Y ( ltt ) 


WRITE 
WRITE 
WRITE 
WRITE 
WRITE 
WRITE 
Y(10) 
Y( 11 ) 

Y ( 12 ) 

Y ( 14 ) 


SAN = 

CAN = 

CXI = 

SXI = 

CXI = 

CAP = 

SAP = 

CAP = 

A ( 202 > = 
A ( 203 ) = 


OUTPUT TAPE 
OUTPUT TAPE 
OUTPUT TAPE 
OUTPUT TAPE 
OUTPUT TAPE 
OUTPUT TAPb 
= Y ( 1 0 ) /RAD 
= Y ( 1 1 ) /RAU 
= Y ( 1 2 ) /RAL) 

= Y ( 14 ) /RAU 
CTHb = THb/RAO 
STHE = SINF(CTHb) 

CTHE = COSE { CTHE ) 

CAN = AN/RAD 

SINF(CAN) 

COSF { CAN ) 

XI /RAD 
SINFKXI ) 

COSE t CXI ) 

( T Hb+AP ) /RAD 
SINF(CAP) 

COSF (CAP) 

SXI*SAP 

SURTF ( 1 • -A ( 202 ) **2 ) 


3 ♦ 1 * AA 9 ECC * X I 
3 * 1 * T Fib * AP 9 AN 
3»1*Y(7)*Y(8)»Y{9) 

3 * 1 ? Y ( 1 0 ) > Y ( 1 1 ) » Y ( 1 2 ) 

3 1 1 9 Y ( 13) » Y ( 14) 

391 9Y ( 15 ) 9 Y<16 ) 9 Y< 17) 9Y< 18 ) 


A ( 204 ) = (SXI*CAP >/A(203) 

A ( 203 ) = CX I / A < 203 ) 

A( 200) = (CAN*CXI*SAP+SAN*CAP ) /A (203) 
A ( 20 1 ) = (CAN*CAP-SAN*CXI*SAP)/A<203) 
P = AA* ( 1 •— ECC**2 ) 

Y ( 4 ) = P/(l.+ECC*CTHt) 

Y ( 5 ) = ATANFt A( 200) *A(201 ) ) 

Y ( 6 ) = ATANF(A(202) *A(203) ) 

Y(l) = SORT F ( ZMU/P ) *ECC*ST HE 
Y ( 3 ) = SORT F ( ZMU* Y ( 4 ) * ( 1 , +ECC*CT HE ) ) 

Y ( 2 ) = Y ( 3 ) *CX I 
Y ( 3 ) = Y(3)*A(204) 

WRITE OUTPUT TAPE 3t2 
CALL RM 14 ) 

END 


SUbROUT INb RK ( N ) 
DIMENSION A(300) »Y( 100) 
COMMON A 9 Y 9 RAD *ZMUfZMU3 
DIMENSION YO ( 50 ) 9 Q ( 50 ) 
INTEGRATION INIAUZATION 
NT = N + 2 
STP = 0. 

DO 1 J= 1 9 N 
1 YU(J) = YU) 

TO o T 
CALL DbRIV 
CALL OUT ( 1 ) 

TPR1NT = Y ( N+l ) + Y ( N+3 ) 
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GO TO 3 

2 00 3 J*1*N 

3 YO(J) = Yt J) 

TO - Y ( N + l ) 

C 1ST STEP Of INTEGRATION 

4 CALL DERI V 

5 DO 6 J= 1 * N 

Y ( J ) = Y0(J)+«5*Y(NT ) *A (J ) 

6 O ( J ) = A ( J ) 

C 2ND STEP Of INTEGRATION 

Y ( N+l ) = Y(N+1)+.3*Y(NT) 

CALL DERIV 
DO 7 J= 1 *N 

YIJ> = Y0tJ)+*3*Y (NT ) *A ( J ) 

7 Q ( J ) = Q ( J ) +2 • * A ( J ) 

C 3RD STEP OF INTEGRATION 

CALL DERIV 
DO 8 J = 1 * N 

Y ( J ) = YO ( J ) + Y ( NT )*A( J) 

8 Q ( J ) = Q ( J ) +2 • *A ( J ) 

C 4TH STEP Of INTEGRATION 

Y ( N+l ) = Y ( N+ 1 ) +• 3* Y ( NT ) 

CALL DERIV 
DO 9 J= 1 * N 

9 Y ( J ) = Y0(J)+(Y(NT)*<U< J)+A( J) ) )/6* 

C CHECK POSSIBLE DAMPER ROD COLLISION 

CALL AUTO( f LAG*STP ) 

If (FLAG) 12*10*12 

10 Y ( N + 1 ) = TO 
DO 11 J = 1 , N 

11 Y ( J ) = YO ( J ) 

GO TO 4 

12 TINT = Y(N+l)-Y(N+4) 

C CHECK TERMINATION 

If ( (TINT/Y(N+4) )+l.E-6) 14,13*13 
C PINAL PRINT BEFORE TERMINATION 

13 CALL OUT ( 1 ) 

GO TO 18 

14 IF ( T I N1 + Y ( NT ) ) 16*16*13 

C GET READY TO TERMINATE AFTER NEXT STEP 

15 Y ( NT ) = -TINT 

C PRINT AND/OR PROCEED 

16 IF ( Y(N+1 J-TPRINT ) 2,17*17 

17 CALL OUT ( 0 ) 

TPRINT = TPR I NT +Y ( N + 3 ) 

GO TO 2 

18 RETURN 
END 


SUBROUTINE DERIV 
DIMENSION A(500) >Y( 100) 

COMMON A , Y , RAD * ZMU * 2MU3 
C TRIGONOMETRIC FUNCTIONS (TABLE D8 ) 

A ( 200 ) = S I NF ( Y ( 5 ) ) 

A ( 201 ) = COSF ( Y ( 3 I ) 

A ( 202 ) = S I NF ( Y ( 6 ) ) 

A l 203) = COSF (Y (6)) 



Table DIO (Continued) 

Program Listing. 

VL = Y(2)/(Y(4)*A(203) ) 

VP = Y(3)/Y(4) 

A ( 205 ) = ATANF ( VP * VL ) 

A ( 204 ) = S I NF ( A < 205 ) ) 

A ( 205 ) = COSF ( A ( 205 ) ) 

A 1 206 ) = SI NF I Y ( 10 ) ) 

A ( 207 ) = COSF ( Y< 10) ) 

A ( 206 ) = SINFIYI 111) 

A ( 209 ) = COSF(Ylll) ) 

A ( 2 1 0 ) = S I NF ( Y ( 12 ) ) 

A ( 21 1 ) = CUSHY (12) ) 

A( 212 ) = SINF ( Y< 14 ) ) 

A ( 2 1 3 ) = COSF ( Y ( 14 ) ) 

A ( 2 14 ) = A( 207 ) *A ( 209 ) 

A ( 2 1 5 ) = -A(207)*A(208)*A(210)-A(206)*A(2ll) 

A(216) = -A(207 >*A(208)*A(211 ) +A l 206 > *A ( 2 10 ) 

A ( 21 7 ) = At 206 ) *A ( 209 ) 

A ( 2 18 ) = -Al206)*A<208)*A<210)+A<207)*A(2ll> 

A ( 2 19 ) = -A(206)*A(208)*A(211)-A(207)*A(2l0) 

A ( 220 ) = A( 209 ) *A( 210 ) 

A ( 22 1 ) = A ( 209 ) *A ( 211) 

A ( 222 ) = A ( 2 13 ) *A ( 65 ) 

A ( 22 3 ) = A(213)*A(64) 

A ( 224 ) = -A( 212 )*A( 65 ) 

A ( 225 ) = -A(212)*A(64) 

R2 = Y ( 4 ) **2 

A ( 100 ) THRU A( 102 ) ARE PROVIDED TO ACCOMMODATE FORCES 
WHICH THE USER MIGHT DLSIR1 TO ADD TO THF PROGRAM 
A ( 100 ) = 0. 

A ( 101 ) = 0. 

A ( 102 ) = 0. 

TRANSLATIONAL DERIVATIVES (TABLE D4) 

All) = A(100)/A(52)+(VL**2+VP**2)/Y(4 )~ZmU/R2 
A ( 2 ) = (A(101)*Y(4)*A(203) ) / A ( 52 ) 

A(3) = (All02)*Y(4) )/A(52)-(VL**2)*(A(202)2A(203) ) 

A(4) = Y(l) 

A ( 5 ) = VL/( Y(4)*A<203) ) 

A ( 6 ) = YI31/R2 

TRIGONOMETRIC FUNCTIONS REQUIRED FOR GRAVITATIONAL MOMENTS 
A ( 226 ) = A( 2 14 ) *A (215) 

A ( 22 7 ) = A( 214 ) *A ( 216 ) 

At 228 ) = A( 2 15 ) *A ( 216 ) 

A ( 229 ) = A ( 2 14 ) **2 
A ( 230 ) = A( 2 1 5 ) **2 
A ( 23 1 ) = A l 2 16 ) **2 
R3 = ZMU3/ ( R2*Y ( 4 ) ) 

ANY ADDITIONAL MOMENTS THE USER MIGHT DESIRE TO ADD 
SHOULD BE INCLUDED IN AI103) THRU A(109) 

A ( 103 ) = R3*(A(228)*A(56)+A(227)*A(59)-A<226>*A(60> 

3+<A(231 )-A(230) )*A(61)) 

A ( 104 ) = R3*< A< 227 ) *A ( 57 ) -A ( 228 )*A ( 59 )+ ( A< 229 )-Al 231 ) ) *A( 60 ) 
4+A 1226) *A (61)) 

A ( 105 ) = R3*(A(226)*A(58)+1A(230)-A(229))*A(59)+A1228)*A(60) 

5-A(227)*A(61 ) ) 

A ( 108 ) = R3*A( 62 )*( A( 214)*A( 212 )+A ( 2 15 ) *A< 222 ) +A( 216 ) *A( 223 ) ) 

A ( 107 ) = A(108)*(A(214>*A(213)+A(215)*A(224)+A(216)*A(225I ) 

A ( 108 ) = -A< 108)*(-A(215)*A(64)+A(216)*A(65) ) 

A ( 109 ) = A(66)*Y ( 13 )+A( 67 )*Y( 14) 

A(110) THRU A ( 1 12 ) ~ ANGULAR MOMENTUM 

A ( 1 10 ) = A ( 53 ) *Y ( 7 ) +A (59)*Y(8)+A(60)*Yt> ; 



n n 
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C 

C 

C 

C 

C 

C 


A ( 1 1 1 ) = A(59)*Y(7)+A(54)*Y(8)+A(61)*Y(9) 

A ( 1 12 ) = A(60)*Y17>+A(61)*Y(8)+A(55)*Y(9) 

A < 1 13 ) THRU A ( 115) - INERTIAL BODY RATES 

A 1 1 1 3 ) = Y(7)*A(212)+Y(8)*A( 222)+YI9)*A(223) 

AC 1 1 A ) = -YC8)*A<64)+Y(9)*AC65)-Y< 13) 

A( 115 J = Y(7)*A>2 ,;, ' J-V 3)*A(224) + Y(9)*A(225) 

A ( 1 16 ) - COUPLING TERM btTWEEN ROD AND MAIN BODY 
M(116> = A ( 62 ) *A C 1 13 ) * ( A ( 1 14 ) - Y ( 1 3 ) ) 

AC 117 ) THRU A ( 119) “ MOMENTS GROUPED TOGETHER IN 

THE ANGULAR RATE EQUATIONS (SEE EQ. 12) 

A ( 1 17 ) = A(103)+A(106)*A(212)+A(108)*A(2l3)“Y(8)*A(112) 

7+Y(9)*A(lll>-A( 116)*A(213) 

A ( 1 1 8 ) = AC 1041+AC 106)*A<222 )+A( 109)*A(64>+A< 108)*AC 224) 

8- Y C 9 ) *A ( 1 10 ) +Y C 7 ) *A C 1 12 ) — A ( 1 16 ) # A ( 224 ) 

AC 119) = A<105)+A<106)*A(223)-A<109>*AC65>+A<108)*A<225) 

9- Y ( 7 ) *A ( 111 ) + Y ( 8 ) *A ( 110 ) —At 1 16 ) *A ( 225 ) 

A C 120 ) THRU AC 122) - ORBITAL OR TRANSLATIONAL RATES 

A( 120 ) = ( VL*A( 102 )-VP*A( 101 ) )/( < VL**2+VP**2 ) *A ( 52 ) ) 

A ( 12 1 ) = 0. 

A ( 122 ) = A(5)*A(203)*A(205)+A(6)*A<204) 

A ( 123 ) THRU AC 125) - ANGULAR RATES RELATIVE TO ORBIT 

AC 123 ) = YC7J-AC 120 ) *A ( 2 14 ) -AC 122 > *A ( 208 ) 

A ( 124 ) = Y(8)-A(120)*A(215)-A(122)*Al220) 

A C 125 ) = Y(9)-A(120)*A<216>-A(122)*A<221) 

A ( 126 ) AND AC 127) - DAMPER ROD JR IGONOMeTRIC FUNCTIONS 

A ( 126 ) = A ( 2 12 ) **2 
A ( 127 ) = AC 212 ) *A ( 21 3 ) 

AC 128 ) THRU AC 134) - ELEMENTS AND DETERMINANT 

OF INVERSE INERTIA MATRIX 
A ( 128 ) = AC75 )+A(81 )*A( 126) 

A ( 129 ) = AC76I+AC82 )*A( 127)+A(83)*A( 126) 

A 1 1 30 ) = AC77)+A(84)*A( 127 > +A ( 85 ) *A < 126) 

A C 1 34 ) = A 1 62 ) *A (213) 

A ( 134 ) = AC128)*(A(53)+A(134)*A(213) ) +A ( i29)*(A(59) 

4 + A ( 134 ) *A( 224 ) )+A( 130)*( AC 60 )+A( 134 )*AC22& ) 1 
A ( 128 ) = A(128)/A( 134) 

A ( 129 ) = A ( 129 ) /A ( 134 ) 

AC 130 ) = AC 130) /A( 134) 

AC 131) = ( A ( 78 )+A(86)*A< 12 7 ) +A ( 8 7 ) *A < 126) ) /AC 134) 

A ( 1 32 ) = (A(79)+A(88)*A( 126 ) +A 1 89 ) *A l 127) ) /AC 134) 

AC 133) = ( A l 80 ) +A ( 90 ) *A 1 127 ) +A ( 91 ) *A ( 126 ) ) /AC 1 34 ) 

ROTATIONAL DERIVATIVES (TABLE 4) 

A ( 7 ) = AC 117 ) * A ( 128 )+A( 118)*A< 129) +AC 119)*Al 130 ) 

AC8) = A ( 1 1 7 ) *A ( 129)+A( 1 18 ) *A ( 131 )+A< 1 1 9 ) * AC 132 ) 

A ( 9 ) = AC 117 ) * A ( 130 )+A( 118)*A( 132 )+AC 119)*A( 133) 

A ( 10 ) = ( AC 124 )*A( 210 )+A( 125 )*A( 211 ) ) /AC209 ) 

A ( 1 1 ) = -A(124)#A(211) +A < 125 l*A < 2 1 0 ) 

A ( 12 ) = AC 123)-A( 10)*At 208) 

A ( 1 3 ) = -AC 8 )*A (64 ) +AC 9 )*A( 65 ) —A ( 1 13 ) *A ( 1 15 )- ( AC 107 ) 

3+AC109) ) / A ( 62 ) 

A ( 14 ) = Y ( 1 3 ) 

RETURN 

END 


SUBROUTINE AUTOC FLAG»STP ) 
DIMENSION AC500) »Yl 100) 
COMMON A»Y»RAD.ZMU»ZMU3 


51 
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RH2 = ABSF( Y( 14) ) 

IF (STP) 1*1,3 

1 IF (RH2-1.) 7*7*2 

SAVE OLD VALUE OF INTEGRATION TIME INTERVAL 

2 DT = Y< 16) 

Y ( 16 ) = Y ( 16 ) * ( ( 1 .-RH1 ) / (RH2-RH1 ) ) 

FLAG = 0. 

STP = 2. 

GO TO 8 

3 STP = STP-1 • 

IF (STP) 4.4,5 

RESTORE OLD VALUE OF INTEGRATION TIME INTERVAL 

4 Y ( 16 ) = DT 
GO TO 7 

5 Y ( 16 ) = DT-Y ( 16 ) 

IF (Y( 13)*Y( 14) ) 7,7,6 
MOMENTUM TRANSFER (TABLE D7 ) 

6 A ( 1 35 ) = A( 110 )+A(62 ) *A( 115 )*A( 213 ) 

A ( 1 36 ) = A(lll )-A<62 )*( A( 1 14 ) *A ( 64 ) -A ( 1 15 > *A ( 224 ) -Y ( 13)*A<64> ) 

A ( 1 3 7 ) = A( 112 >+A(62 )*( Al 114 )*A(65 )+Al 1 15 )*A ( 225)-Y ( 13 )*Al 65 ) ) 

A ( 1 38 ) = A1128)*A(134)+A(97)*A<126)+A<92) 

A ( 1 39 ) = A(129)*A(134>+A(97)*A(222)*A1212>+A(93)*A(65) 

A ( 140 ) = A(130)*A(134)+A(97)*A(223)*A(212) + A(93)*A(64) 

A ( 144 ) = A ( 62 ) *A ( 2 1 3 ) 

A ( 144 ) = A(138)*(A(53)+A( 144 ) *A ( 2 1 3 ) ) +A ( 1 39 ) * ( A I 59 ) 

4+ A ( 144 ) *A t 224 ) )+A( 140 ) * l A 1 60 ) +A l 144>*A(22&> ) 

4+ A ( 140 ) *( A ( 60 ) +A ( 144 ) *A ( 225 ) ) 

A ( 1 38 ) = A ( 1 38 ) / A ( 144 ) 

A ( 1 39 ) = A ( 1 39 ) / A ( 144 ) 

A ( 140 ) = A ( 140 ) /A ( 144 ) 

A ( 141 ) = (A( 131 )*Al 134)+A(y4)+A{97)*A(222) **2 ) / A ( 144 ) 

A ( 142 ) = ( A< 132)*A< 134)+A(95 1+AI97 )*A(222 >*A< 223) ) /A< 144 ) 

A ( 143 ) = (A(133)*A(134)+A(96)+A(97)*A(223)**2>/A(144) 

YI7) = A( 1 35 ) *A ( 138) +A( 1 36 ) *A ( 1 39 ) +A ( 13 7 ) *A ( 140) 

Y(6) = A( 135 )*A( 139 )+A( 136 )*A< 141 ) +A ( 137 ) *A< 142 ) 

Y ( 9 ) = A( 135)*A( 140 )+A( 1 36 ) *A ( 142 ) +A ( 137)*A(143) 

20 FORMATI8X13H HIT STOP AT Ell. 4, 8H RAD/SEO 
WRITE OUTPUT TAPE 3,20*Y(13) 

Y( 13) = -Y (13) 

7 RH1 = KH2 
FLAG = 1. 

8 RETURN 
END 


SUBROUTINE OUT (Id 
DIMENSION A ( 500 ) ,Y< 100) 
COMMON A, Y »RAD ,ZMU *ZMU3 
THR = Y ( 1 5 ) / 3600. 

PITCH = Y ( 10 ) *RAD 
ROLL = Y ( 1 1 ) *RAD 
YAW = Y ( 12 ) *RAD 
INT = YAW/360. 

XINT = INT 

YAW = Y AW-X I NT*360 « 

RH02 = Y(14)*RAD 
TEST OUTPUT MODE 
IF (X) 1*1*2 
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OUTPUT MODE (0) (TABLE D3 ( 

1 WRITE OUTPUT TAPE 3 , 20 * THR »P i TCH.ROLL , Y Aw »RH02 
2u PORMAT ( 5F 12 ,4 ) 

00 TO 9 

2 CXI = A ( 203 ) *A ( 205 ) 

SX I = SQRTF ( l.-CXI**2 ) 

XI = ATANF(SXI *CXI )*RAD 
IF (XI) 3.4,3 

3 SAN = A(200)*A(204)-A(201)*A(202)*A<205) 

CAN = A(201 )*A(204)+A(200)*A(202)*A(205) 

AN = AT ANF ( SAN « CAN ) *RAD 

CAP = A ( 203 ) *A ( 204 ) 

AP = ATANF( A(202 ) ,CAP)*RAD 
GO TO 5 

4 AN = 0. 

AP = ATANF(A(200) ,A(201 ) )*RAD 

5 VL = Y(2)/(Y(4)*A(203) ) 

VP = Y ( 3 ) / Y ( 4 ) 

VT = VL**2+VP**2 
V = Y ( 1 ) **2+VT 
Q = (Y(4)*V)/ZMU 
CL2 = VT/V 
AA = Y(4)/(2.-Q) 

ECC = SQRTF ( l.-Q* ( 2 .-Q ) *CL2 ) 

IF (ECC-l.E-5) 7,7,6 

6 CTHE = SQRTF(VT)-SQRTF(ZMU/(AA*(1.-ECC**2) ) ) 
THE = ATANF ( Y( 1 ) »CTHE )*RA0 

AP = AP-THE 
GO TO 8 

7 THE = AP 
AP = 0. 

ECC = 0. 

8 V = SQRTF (V) 


OUTPUT 

MODE ( 

1 ) 

(TABLE D3 ) 


WRITE 

OUTPUT 

TAPE 

3,23 


WRITE 

OUTPUT 

TAPE 

3,21, THR , Y ( 15 ) 


WRITE 

OUTPUT 

TAPE 

3 , 22 , P I TCH , Y ( 7 ) , AA, 

THE 

WRITE 

UUTPUT 

TAPE 

3,22, ROLL , Y ( 8 ) .ECC, 

AP 

WRITE 

OUTPUT 

TAPE 

3,22 , YAW , Y ( 9 ) ,XI »AN 


WRITE 

OUTPUT 

TAPE 

3 , 21 , Y ( 4 ) , V,RH02 ,Y ( 

13) 

WRITE 

OUTPUT 

TAPE 

3,23 



21 FORMAT ( 4E18 .8 ) 

22 FORMAT (E20,8,3E18.8) 

23 FORMAT (1H0) 

9 RETURN 

END 



